rm(list = ls())

wd <- ".../Replication/"
setwd(wd)

# Load/install packages --
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  foreign, 
  ggplot2, 
  estimatr,
  texreg,
  xtable,
  fastDummies,
  sandwich,
  dplyr,
  janitor, 
  gridExtra,
  gsheet,
  zoo,
  interflex,
  lubridate,
  tidyverse,
  stringi,
  readxl,
  ri2,
  modelsummary,
  ggpubr
)

options(scipen=999)

# Note: given the random component inherent to the randomization tests performed
# by ri(), you may observe small differences in the p-values of the randomization
# tests performed using that function.

# Load data ----
d <- read.csv("Data/baseline database.csv")
winners <- subset(d, d$group=="Liberado sorteo")
censo <- read.csv("Data/census_nonwhites.csv")

# remove from the census individuals with possible matches in treated group
censo$name1 = paste(censo$nombre,censo$apellido)
censo$name2 = paste(censo$nombre,censo$apellido_jefehogar)
censo$name_search = ifelse(is.na(censo$apellido)==TRUE,censo$name2,censo$name1)
censo$name_search = ifelse(is.na(censo$apellido)==TRUE & is.na(censo$apellido_jefehogar)==TRUE,NA,censo$name_search)
censo = subset(censo, is.na(censo$name_search)==FALSE)

possiblematches = NA

for(i in 1:nrow(winners)){
  list1 = agrep(winners$full_name[i], censo$name_search, max.distance =0.10, value = TRUE)
  possiblematches = c(possiblematches,list1)
}

censo2 <- subset(censo, !(censo$name_search %in% possiblematches))
enslaved <- subset(censo2, censo2$slave==1)
enslaved$lottery_winner <- 0
enslaved$imprisoned_dummy <- enslaved$imprisoned>0

win <- cbind.data.frame(winners$full_name,
                        winners$imprisoned_dummy,
                        winners$lottery_winner)


# create subsets of the census
all <- cbind.data.frame(enslaved$nombre,
                        enslaved$imprisoned_dummy,
                        enslaved$lottery_winner)

all_15_20 <- enslaved %>% filter(age>=15 & age<=20) %>% 
  select(nombre,imprisoned_dummy,lottery_winner)
all_15_30 <- enslaved %>% filter(age>=15 & age<=30) %>% 
  select(nombre,imprisoned_dummy,lottery_winner)
all_15_40 <- enslaved %>% filter(age>=15 & age<=40) %>% 
  select(nombre,imprisoned_dummy,lottery_winner)

colnames(win) <- c("name","imprisoned_dummy","lottery_winner")
colnames(all) <- c("name","imprisoned_dummy","lottery_winner")
colnames(all_15_20) <- c("name","imprisoned_dummy","lottery_winner")
colnames(all_15_30) <- c("name","imprisoned_dummy","lottery_winner")
colnames(all_15_40) <- c("name","imprisoned_dummy","lottery_winner")

# final data
d1 <- rbind.data.frame(win,all)
d2 <- rbind.data.frame(win,all_15_20)
d3 <- rbind.data.frame(win,all_15_30)
d4 <- rbind.data.frame(win,all_15_40)

# compute exact p-values
declare.d1 <- declare_ra(N = nrow(d1), m=sum(d1$lottery_winner)) 
declare.d2 <- declare_ra(N = nrow(d2), m=sum(d2$lottery_winner)) 
declare.d3 <- declare_ra(N = nrow(d3), m=sum(d3$lottery_winner)) 
declare.d4 <- declare_ra(N = nrow(d4), m=sum(d4$lottery_winner)) 

#---- Tables ----

# estimate regressions
ITT <- list()
ITT[['all']] <- lm_robust(imprisoned_dummy ~ lottery_winner, data = d1)
ITT[['all_15_40']] <- lm_robust(imprisoned_dummy ~ lottery_winner, data = d4)
ITT[['all_15_30']] <- lm_robust(imprisoned_dummy ~ lottery_winner, data = d3)
ITT[['all_15_20']] <- lm_robust(imprisoned_dummy ~ lottery_winner, data = d2)

ri <- list()
ri[['all']] <- summary(conduct_ri(
  formula = imprisoned_dummy ~ lottery_winner,
  declaration = declare.d1,
  sharp_hypothesis = 0,
  assignment = "lottery_winner",
  data = d1,
  sims = 10000
))


ri[['all_15_40']] <- summary(conduct_ri(
  formula = imprisoned_dummy ~ lottery_winner,
  declaration = declare.d4,
  sharp_hypothesis = 0,
  assignment = "lottery_winner",
  data = d4,
  sims = 10000
))

ri[['all_15_30']] <- summary(conduct_ri(
  formula = imprisoned_dummy ~ lottery_winner,
  declaration = declare.d3,
  sharp_hypothesis = 0,
  assignment = "lottery_winner",
  data = d3,
  sims = 10000
))


ri[['all_15_20']] <- summary(conduct_ri(
  formula = imprisoned_dummy ~ lottery_winner,
  declaration = declare.d2,
  sharp_hypothesis = 0,
  assignment = "lottery_winner",
  data = d2,
  sims = 10000
))


Table_3<-modelsummary(ITT, stars = T, output = "dataframe")
Table_3$part<-NULL
Table_3[nrow(Table_3)+1,1]<-'Exact p-value'
Table_3[nrow(Table_3),'all']<-ri[['all']][,"two_tailed_p_value"]
Table_3[nrow(Table_3),'all_15_40']<-ri[['all_15_40']][,"two_tailed_p_value"]
Table_3[nrow(Table_3),'all_15_30']<-ri[['all_15_30']][,"two_tailed_p_value"]
Table_3[nrow(Table_3),'all_15_20']<-ri[['all_15_20']][,"two_tailed_p_value"]
View(Table_3)

print(xtable(Table_3), file = "Output/Table 3.tex")



